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Abstract: The status of regional biodiversity is determined by habitat quality. The effective assessment of 
habitat quality can help balance the relationship between economic development and biodiversity 
conservation. Therefore, this study used the InVEST model to conduct a dynamic evaluation of the 
spatial and temporal changes in habitat quality of the Tarim River Basin in southern Xinjiang Uygur 
Autonomous Region of China by calculating the degradation degree levels for habitat types that were 
caused by threat factors from 1990 to 2018 (represented by four periods of 1990, 2000, 2010 and 2018). 
Specifically, we used spatial autocorrelation analysis and Getis-Ord G; analysis to divide the study area 
into three heterogeneous units in terms of habitat quality: cold spot areas, hot spot areas and random 
areas. Hemeroby index, population density, gross domestic product (GDP), altitude and distance from 
water source (DWS) were then chosen as the main disturbance factors. Linear correlation and spatial 
regression models were subsequently used to analyze the influences of disturbance factors on habitat 
quality. The results demonstrated that the overall level of habitat quality in the TRB was poor, showing a 
continuous degradation state. The intensity of the negative correlation between habitat quality and 
Hemeroby index was proven to be strongest in cold spot areas, hot spot areas and random areas. The 
spatial lag model (SLM) was better suited to spatial regression analysis due to the spatial dependence of 
habitat quality and disturbance factors in heterogeneous units. By analyzing the model, Hemeroby index 
was found to have the greatest impact on habitat quality in the studied four periods (1990, 2000, 2010 and 
2018). The research results have potential guiding significance for the formulation of reasonable 
management policies in the TRB as well as other river basins in arid areas. 
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1 Introduction 


Biodiversity is the combination of an organism, complex environment and various related 
ecological processes and is also the material basis for human survival (Ma et al., 1998). As human 
activities increase and the economy and society develop rapidly, human disturbance to ecosystems 
has affected the quality of biological habitats and caused a reduction in biodiversity. According to 
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statistic data (2011-2020), over 1.5x10* species have disappeared, and the trend of biodiversity 
loss has not been effectively curbed; biodiversity is anticipated to decrease further in the future 
(Isbell et al., 2013; Terrado et al., 2016a; Dai et al., 2018). Habitat quality is the ability of an 
ecosystem to provide suitable living conditions for the sustainable development of individuals and 
populations within a certain space-time range, which reflects regional biodiversity status to a 
certain extent (Fellman et al., 2015; Hillard et al., 2015). Therefore, the exploration of the 
changes in regional habitat quality and analysis of the influencing factors of these changes are of 
great significance for the protection and ecological management of regional biodiversity. 

Habitat quality is regarded as an important representation of regional biodiversity and 
ecosystem services, and the assessment, simulation and prediction of its trend and status are an 
effective means for studying regional ecosystem services (Wu et al., 2021). Currently, the main 
methods used for habitat quality assessment include the traditional biodiversity and habitat survey 
method (Mark et al., 2008; Miller et al., 2010), ecological index evaluation method (Maes et al., 
2012; Coates et al., 2016) and ecological assessment model (Costa et al., 2010; Terrado et al., 
2016a, b). Due to a lack of long-term and continuous species detection data in field biodiversity 
surveys, biodiversity and habitat survey methods cannot be used to evaluate the temporal and 
spatial dynamic changes in biodiversity (Balasooriya et al., 2008; Sun et al., 2010). For 
calculation, the ecological index evaluation method uses the habitat quality evaluation index 
system and evaluation criteria, such as the biological abundance index and vegetation coverage, 
which have certain defects in terms of the evaluation of dynamic changes in habitat quality and its 
spatial agglomeration state. With the wider application of 3S (geographic information system 
(GIS), global positioning system (GPS), and remote sensing (RS)) technology and mathematical 
models, ecological assessment models have become a powerful tool for the quantitative, visual 
and fine monitoring and assessment of temporal and spatial changes in habitat quality (Leh et al., 
2013). Currently, the main ecological assessment models include the Habitat Suitability Index 
(HSI) model (Liu et al., 2006), Social Value for Ecosystem Services (SolVES) model (Sherrouse 
et al., 2014) and Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) model 
(Sharp et al., 2016; Aneseyee et al., 2020). Of these, the InVEST model has the advantages of 
having relatively few operation parameters, easily accessing to basic data and providing 
quantifiable assessment results and spatial visualization; the model is widely used for studying 
habitat quality, particularly in large-scale basins or regions. For example, Chu et al. (2018) 
combined an InVEST model and Cellular Automaton (CA)-Markov method to study the evolution 
of landscape patterns and habitat quality in the Hubei section of the Three Gorges Reservoir Area 
in China. Their results demonstrated that habitat degradation and landscape pattern change are 
causing the continuous decline of biodiversity in the region. Aneseyee et al. (2020) used the 
InVEST model to analyze the temporal and spatial changes in habitat quality of different land use 
types in the Winike Basin (in the Omo-Gibe Basin, Southwest Ethiopia) and discussed the 
influence of population density, altitude, land use intensity and other factors on habitat quality. 

Changes in habitat quality are a result of many factors, whereby variations in land use 
(landscape pattern) caused by anthropogenic factors are an important indicator that threatens 
habitat quality (Zheng et al., 2018; Zhang et al., 2020). In addition, the distribution of habitat 
quality is also affected by natural factors and other anthropogenic factors, including topography, 
terrain, gross domestic product (GDP) and population density. Most previous studies have 
focused on the impact of land use and spatial pattern change on regional habitat quality (e.g., 
Zheng et al., 2018). Researchers have generally directly associated land use change with habitat 
quality models when studying the impact of land use change on habitat quality (Chu et al., 2018; 
Wang et al., 2019; Wang et al., 2020). Studies have also focused on analyzing the correlation 
between habitat quality and natural and anthropogenic factors, including altitude, slope, GDP and 
intensity of land use, as a means of quantifying the impact of different disturbance factors on 
habitat quality (Sun et al., 2019). However, relatively few studies have been conducted on the 
relationship between the comprehensive indicators of human disturbance and habitat quality. 


HE Bing et al.: Assessment of river basin habitat quality and its relationship... 


Furthermore, the heterogeneity of habitat quality has been largely ignored, and the relationship 
between habitat quality and human disturbance can be distorted (Zhang et al., 2018). Previous 
studies (Hu et al., 2015; Han et al., 2020) have demonstrated that the Getis-Ord G; spatial 
statistical method can be used for describing and visualizing the spatial distribution of related 
elements and can effectively identify heterogeneous units that are related to clusters. Therefore, 
conducting research on the spatial heterogeneity of habitat quality and the impact of different 
interference factors on habitat quality can help to reveal the evolution of habitat quality when it is 
influenced by human activities and can scientifically guide the protection of the ecological 
environment and the sustainable use of land resources. 

The Tarim River Basin (TRB) is located in southern Xinjiang Uygur Autonomous Region, 
China. It is a typical inland arid area and also an important area in terms of landscape biodiversity 
protection. The region has an abundance of natural resources and plays a crucial role in strategic 
development, but it has a fragile ecological environment (Xue et al., 2019). The geographic 
components of vegetation found in the TRB include desert vegetation flora. Plant species are 
relatively poor, mainly including Populus euphratica, Tamarix chinensis, Glycyrrhiza uralensis 
and Phragmites communis. The desert riparian forest with Populus euphratica and Tamarix sp. as 
constructive and dominant species is the main barrier for the desertification of oases (Yu et al., 
2017). In recent decades, as a result of increased human activities, social activities and economic 
activities, the unreasonable development and utilization of water and soil resources in the TRB 
have caused several ecological and environmental problems (vegetation degradation, lake 
drought, desertification, massive loss of biodiversity and habitat degradation) (Tan et al., 2011; 
Ling et al., 2018). Following the implementation of the Belt and Road initiative, as the core area 
of the Silk Road Economic Belt (Chen et al., 2016), the habitat quality of the TRB is of great 
significance in terms of social and economic development, ecosystem balance and biodiversity 
protection. However, most existing researches on the TRB have focused on climate change, 
optimization of the allocation of water resources, land use change and ecological water transfer 
effects, whereas comparatively little research was focused on biodiversity from an ecosystem 
service perspective. 

The main objectives of this study are as follows: (1) evaluating habitat quality between 1990 
and 2018 and analyzing the temporal and spatial changes of habitat quality by determining habitat 
suitability and the threat to biodiversity caused by different land use types in the TRB; (2) 
identifying and dividing the spatial heterogeneity units of habitat quality; (3) clarifying the 
differences of influencing factors on habitat quality in different heterogeneous units using linear 
regression method; and (4) clarifying the spatial dependence of habitat quality on disturbance 
factors in different units using a spatial regression model. It is anticipated that the results of this 
study will provide a reference for the rational planning and utilization of land resources and the 
coordinated development of the ecological environment in the TRB. 


2 Materials and methods 


2.1 Study area 


The TRB is located in the Tarim Basin (34°55'-43°08’N, 73°10’—94°05’E) in southern Xinjiang 
Uygur Autonomous Region of China and covers an area of 1.02x10° km’, which includes 9 water 
systems and 114 rivers (Yu et al., 2016). The TRB is located at a significant distance from the 
ocean and is surrounded by mountains on every side. It is a closed inland hydrological region 
with a relatively independent water cycle and water balance (Fig. 1) (Xue et al., 2017). The Tarim 
River is China's largest inland river with a mainstream that is 1321 km long and does not produce 
any flow. Historically, the water from nine water systems in the TRB flowed into the mainstream 
of the Tarim River, but due to a change in the natural environment and the development of 
modern oasis agriculture, only three upstream rivers (Aksu River, Yarkant River and Hotan River) 
now have a natural hydraulic connection with it (Chen et al., 2013; Xu et al., 2013). The TRB 
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encompasses 5 states, 42 counties (cities), 4 production and construction companies and 55 
regiments. According to statistical data, in 2018, the total population of the TRB was 13.14x10°, 
the crop planting area was 2.23x10° hm’, the gross national product was 411.77x10° CNY and the 
total water resources were 42.0410? m°. 


JERE 


Legend Elevation (m) 

— River ; High: 8439 

© Basin boundary — Low: 41 
Fig. 1 Overview of the Tarim River Basin (TRB). ARB, Aksu River Basin; YRB, Yarkant River Basin; HRB, 
Hotan River Basin; CRB, Cherchen River Basin; WRB, Weigan-Kuche River Basin; KRB1, Kaxgar River Basin; 
KRB2, Keriya River Basin; KRB3, Kaidu-Konqi River Basin. 
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2.2 Data sources and preprocessing 


In this study, we used five different data sources to calculate and analyze the temporal and spatial 
variations in habitat quality and spatial heterogeneity with disturbance factors. Specifically, 
Landsat series RS images with a resolution of 30 mx30 m from 1990, 2000, 2010 and 2018 in the 
TRB were chosen as the basic data sources, and the datasets were provided by the Data Centre for 
Resources and Environmental Sciences at the Chinese Academy of Sciences (RESDC) 
(http://www.resdc.cn). RS images with no clouds or few clouds between June and September 
were selected. Using ENVI 5.3 software, radiation correction, geometric correction, clipping and 
other preprocessing operations were conducted. The land use information was extracted and 
assigned using the human-computer interaction method, and the interpretation results were 
verified and corrected using GPS technology. Finally, the land use thematic map was obtained 
with an accuracy of 97.15% and unified classification. Land use types were categorized as 
cultivated land, forestland, grassland, water body, construction land and unused land using 
ArcGIS 10.5 software (Fig. 2) (Sun et al., 2021). The 90 mx90 m DEM (digital elevation model) 
data from the study area were obtained from the geospatial data cloud (http://www.gscloud.cn/). 
The 1 kmx1 km GDP data (1990, 2000, 2010 and 2018) were obtained from the RESDC. River 
system data (1:1,000,000) were obtained from the RESDC. Population data (1991-2019) were 
obtained from the Xinjiang Statistical Yearbook (Statistics Bureau of Xinjiang Uygur 
Autonomous Region, 1990-2018). 


2.3 Model and methods 
2.3.1 Model 


The InVEST model was developed by the US natural capital project team in 2007, and the habitat 
quality module is the one that has multiple ecosystem services for assessment (Sharp et al., 2016). 
The habitat quality module reflects the impact of human activities on the ecological environment; 
therefore, the greater the disturbance of human activities is, the greater the threat to the habitat, 
and the lower the quality is, the lower the biodiversity. Habitat type and threat source factors must 
be set when the habitat quality module is running (Zhang et al., 2011). Due to the different 
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Fig. 2 Temporal and spatial distributions of land use types in the TRB in 1990 (a), 2000 (b), 2010 (c) and 2018 (d) 


sensitivities of various land use types to threat source factors (Yang et al., 2018), the current 
situation and expert opinions of the TRB were considered in this study, and cultivated land, 
construction land and unused land were chosen as threat source factors. The other land use types 
(forestland, grassland and water body) represent different habitat types. After referring to the user 
manual for the InVEST model and examining previous research results (Zhang et al., 2011; Liu et 
al., 2020; Wang et al., 2021), we then determined the weight of threat source factors, the 
maximum influence distance of source factors on ecological land (Table 1) and the sensitivity of 
ecological land to stress factors (Table 2). The core objective of the habitat quality module is to 
establish the relationship between habitat types and threat source factors or to obtain the 
degradation degree of habitat by calculating the negative impact of threat source factors on the 
habitat and then determining the habitat quality based on the suitability and degradation degree of 
habitat. Before the model was run, all land use data and threat source factors data were rasterized 
using ArcGIS software. The specific calculation equations are as follows (Han et al., 2019): 


r Y, W, i 
D=} > R PTN (1) 
r=l y=1 $w, 
r=l 


r max 


d 
iy =1- í = Jor linear); i, = exp -( 22 Je a Gf exponential), (2) 


Q;=H;|1-| —— ||, (3) 


where j is the habitat type; Dx is the habitat degradation degree of grid x in habitat type j; r is the 
number of threat source factors (r=3 in this study); Y, is the total number of grids of threat source 
factor r; W, is the weight of threat source factor r; ry is the stress value of threat source factor r to 
grid y; Py is the accessibility of grid y; S; is the sensitivity of habitat type j to the threat source 
factors r, and the closer the value is to 1, the more sensitive it is; i-xy is the stress value of grid y; 
dxy is the distance between grid x and grid y ;d, is the influence range of threat source factor r; Qx 
is the habitat quality index of grid x in habitat type j; Hj is the habitat suitability of habitat type j; 
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K is the semi-saturation constant, generally, half of the maximum value of D,j; and z is the 
normalization constant, which is taken as 2.5 by default. 


Table 1 Attribute table of habitat threat source factors 


Threat source factor (R) dimax (km) Weight Spatial decay function 
Cultivated land 0.35 0.60 Linear 
Construction land 8.00 0.40 Exponential 
Unused land 3.00 0.50 Linear 


Note: d,max is the maximum influence range of threat source factor r. 


Table 2 Sensitivity of land use types to each threat source factor 


Sensitivity 
Code Land use type 
Habitat suitability Cultivated land Construction land Unused land 
10 Cultivated land 0.30 0.00 0.90 0.50 
20 Forestland 1.00 0.60 0.80 0.20 
30 Grassland 0.90 0.80 0.50 0.30 
40 Water body 1.00 0.50 0.40 0.50 
50 Construction land 0.00 0.00 0.00 0.00 
60 Unused land 0.10 0.10 0.30 0.00 


To compare and analyze the habitat quality changes in different years of the TRB, we used the 
natural break (jenk) classification in ArcGIS with four grades: low (0.00-0.25), lower 
(0.25—0.50), higher (0.50-0.75) and high (0.75—1.00). 


2.3.2 Methods 


Due to the large area of the TRB, we divided the study area into 113,916 grid units with a spatial 
resolution of 3 kmx3 km using the "create fishing net" function module of ArcGIS 10.5 software, 
to analyze the spatial heterogeneity of habitat quality and its relationship with disturbance factors. 
Subsequent spatial autocorrelation analysis, cold-hot spot analysis, ecological disturbance degree, 
and the relationship between habitat quality and disturbance factors were all based on the created 
grids. The software used for the calculation included ArcGIS 10.5, GeoDa 1.1 and SPSS 24.0. 
Hemeroby index is a comprehensive indicator of human disturbance and is determined on the 

basis of the comprehensive impact human activity frequency and the degree of human activities 
on an ecosystem (Hill et al., 2002; Chen et al., 2010). The calculation equation is as follows: 

> HI,;xS, 

i=l 


f=. 4 
x (4) 


where H is the Hemeroby index of a grid cell; HI; is the Hemeroby coefficient of the i® land use 
type (Table 3); S; is the i land use type in a grid cell; and S is the total number of the grid units. 
Based on the natural break point, we divided the Hemeroby index into three levels: low 
(0.0000—0.4000), medium (0.4000—0.6000) and high (0.6000—1.0000). 


Table 3 Hemeroby coefficient (HI) corresponding to each land use type 


Cultivated land Forestland Grassland Water body Construction land Unused land 


HI 0.7500 0.5500 0.5000 0.3000 0.9500 0.7200 


Spatial autocorrelation analysis is a method that is used for testing the potential dependence of 
spatial variables with a certain degree of regularity in different spatial positions based on classical 
Statistics (Rey, 2001). As an important spatial statistical research field, the combination of global 
spatial autocorrelation analysis and local spatial autocorrelation analysis has been widely used in 
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a variety of spatial problems (Zhang et al., 2011). Spatial heterogeneity is resulted from different 
degrees of spatial autocorrelation (Han et al., 2020). Therefore, this study used the global spatial 
autocorrelation analysis method to test the heterogeneity of habitat quality within the study area. 
Generally, the global Moran's J index was used for measurement, and the calculation equation is 
as follows: 


n> Wii x) 
i= n = = n d (5) 
Wy æy 
i=l j=1 i=l 
where J is the Moran's J index; n is the total number of grids; W; is the spatial weight matrix of 
raster i and raster j (if i is adjacent to j, its spatial weight is 1; otherwise, it is 0); x; is the habitat 
quality values in the i grids; and x is the mean value of habitat quality. 

When conducting spatial autocorrelation analysis, the habitat quality values for the four periods 
were extracted using 3 kmx3 km grids. They were then processed and calculated using GeoDa 
software. 

The method for studying the clustering distribution characteristics of local regions is known as 
hotspot analysis. This is a method that is used for identifying statistically significant high-value 
areas (hot spots) and low-value areas (cold spots) in the spatial distribution of habitat quality (Han 
et al., 2019). In this study, ArcGIS software was used to analyze the habitat quality of the TRB in 
1990, 2000, 2010 and 2018 using the Getis-Ord G; index (Getis et al., 1992). According to the 
calculation results of the Getis-Ord G; index, we divided the study area into cold spot areas, hot 
spot areas and random areas. Compared with the Moran's 7 index, the Getis-Ord G; index is able 
to determine the spatial distribution of habitat quality heterogeneity units. The calculation 
equation is as follows: 


(6) 


where x; is the habitat quality of raster j; and S is the standard deviation of habitat quality. 

Taking the 3 kmx3 km grid element as the calculation scale, the Pearson correlation 
coefficients of habitat quality and disturbance factors in the three habitat quality heterogeneous 
units (cold spot areas, hot spot areas and random areas) were calculated for 1990, 2000, 2010 and 
2018 using the linear regression tool of SPSS software. The three habitat quality heterogeneity 
units were determined using cold-hot spot analysis. 

In this study, the ordinary least square (OLS), spatial lag model (SLM) and spatial error model 
(SEM) proposed by Anselin (1988) were used to perform regression analysis of habitat quality 
and disturbance factors using GeoDa software. The calculation equation for the model is as 
follows (Anselin, 2005): 

y=pW, y+ hx+u+ au = AW, +E, (7) 


where y is the dependent variable, i.e., habitat quality; p is the regression coefficient of spatial lag 
term Wiy; Wı and W2 are the spatial adjacency weight matrices of the dependent variable and 
residual, respectively; # is the regression coefficient of the independent variable; x is the 
independent variable, i.e., the relevant interference factor; u is a random error term; æ is a 
constant; J is the regression coefficient of the spatial residual term; and ¢ is a random error with a 
mean value of 0 and variance of 67. 

When different parameters in Equation 7 are equal to 0, four types of spatial regression models 
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can be formed (Anselin, 2005). When p40 and f=A=0, this is a first-order spatial autoregressive 
model that does not take the influence of the independent variables on the dependent variables 
into account. Therefore, only the OLS, SEM and SLM models were considered for use in this 
study. When p=0 and A=0, the OLS model is defined, and the observed values of the dependent 
and independent variables in space are unaffected by spatial differences. When p40 and /=0, the 
SLM model is defined, and there is a spatial correlation between the dependent variables in space, 
with the observed values of the dependent variables being related to the corresponding 
independent variables and the dependent variables in the adjacent areas. When p=0 and 40, the 
SEM model is defined. In this case, there is no spatial correlation between the dependent 
variables in space, and only the independent variables with spatial correlation are considered. In 
addition, the observation values of the dependent variables in space are related to the 
corresponding independent variables as well as the independent variables and dependent variables 
in adjacent areas. 


3 Results 


3.1 Temporal and spatial variations of habitat quality in the TRB 


Based on the dynamic change in land use types in the TRB between 1990 and 2018, we evaluated 
habitat quality of four periods (1990, 2000, 2010 and 2018) using the habitat quality module of 
the InVEST model. In ArcGIS, habitat quality values were divided into four grades: 
0.0000-0.2500, 0.2500-0.5000, 0.5000—0.7500 and 0.7500-1.0000. The smaller the value was, 
the worse the habitat quality. Figure 3 shows that the overall level of habitat quality for the entire 
TRB was poor (in the range of 0.2500—0.5000), and the average habitat quality continued to 
decrease from 0.3000 in 1990 to 0.2700 in 2018, which indicated that habitat quality in the study 
area was in a State of degradation. Between 1990 and 2018, habitat quality in each sub-basin 
exhibited different degrees of degradation (with the exception of the Hotan River Basin (HRB)). 
The spatial distribution of habitat quality demonstrated that there were differences in habitat 
quality between sub-basins in the study area (Tables 4 and 5). Habitat quality of the northern 
sub-basins (Kaxgar River Basin (KRB1), Aksu River Basin (ARB), Weigan River Basin (WRB) 
and Kaidu-Konqi River Basin (KRB3) was higher than that of the southern sub-basins (HRB, 
Keriya River Basin (KRB2), and Cherchen River Basin (CRB)), and habitat quality of 
high-altitude areas (gradient 3 (3000—4500 m) and gradient 4 (4500-8439 m); Table 5) was higher 


(a) 1990 p (b) 2000 2 À 
(c) 2010 | > (d) 2018 | > 


Legend 
Habitat quality 0 
E 0.0000-0.2500 mı 0.2500—-0.5000 mm 0.5000-0.7500 mm 0.7500—1.0000 


400 800 km 


Fig.3 Spatial distribution of habitat quality in the TRB in 1990 (a), 2000 (b), 2010 (c) and 2018 (d) 
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Table 4 Changes in habitat quality in the sub-basins of the TRB from 1990 to 2018 
Habitat quality 


Sub-basin 
1990 2000 2010 2018 Value change (1990-2018) 
ARB 0.4476 0.4362 0.4292 0.4146 —0.0330 
KRB1 0.4530 0.4552 0.4156 0.4097 —0.0433 
YRB 0.4410 0.4391 0.4282 0.4257 0.0153 
HRB 0.3605 0.3521 0.3729 0.3724 0.0119 
KRB2 0.3597 0.3533 0.2894 0.2896 0.0701 
CRB 0.2878 0.2895 0.2715 0.2713 —0.0165 
KRB3 0.3649 0.3634 0.3327 0.3327 —0.0322 
WRB 0.4268 0.4290 0.3649 0.3584 —0.0684 
Mainstream 0.4968 0.4828 0.4471 0.4459 —0.0509 
Taklimakan Desert 0.1409 0.1365 0.1203 0.1205 —0.0204 
Kumtag Desert 0.1384 0.1371 0.1351 0.1342 —0.0042 
TRB 0.2981 0.2950 0.2766 0.2746 —0.0235 


Note: TRB, Tarim River Basin; ARB, Aksu River Basin; KRB1, Kaxgar River Basin; YRB, Yarkant River Basin; HRB, Hotan River 
Basin; KRB2, Keriya River Basin; CRB, Cherchen River Basin; KRB3, Kaidu-Konqi River Basin; WRB, Weigan-Kuche River Basin. 


Table 5 Habitat quality changes in the TRB at different altitudinal gradients from 1990 to 2018 
Habitat quality 


Altitudinal gradient 


1990 2000 2010 2015 Value change (1990-2018) 
Gradient 1 0.2266 0.2222 0.1991 0.1980 —0.0330 
Gradient 2 0.3377 0.3381 0.3311 0.3317 —0.0433 
Gradient 3 0.5238 0.5229 0.4802 0.4743 —0.0204 
Gradient 4 0.3727 0.3697 0.3795 0.3775 —0.0042 
TRB 0.2981 0.2950 0.2766 0.2746 —0.0235 


Note: Gradient 1, 41-1500 m; Gradient 2, 1500-3000 m; Gradient 3, 3000-4500 m; Gradient 4, 4500-8439 m. The division of elevation 
gradient 1 is based on the fact that most hydrological stations at the outflow of the mountains are 1500 m above sea level, and the 
outflow hydrological station is divided into mountainous and plain areas. In addition, the classification of elevation gradient 4 is based 
on the existence of permanent snow and glaciers in the mountainous areas of above 4500 m. Elevation gradient 2 and elevation gradient 
3 are divided on the basis of the isometric interval. 


than that of low-altitude areas (gradient 1 (41-1500 m) and gradient 2 (1500-3000 m); Table 5). 
The spatial differences of habitat quality between 1990 and 2018 were obvious. With the 
exceptions of the Taklimakan Desert and Kumtag Desert, the lowest habitat quality was found to 
be mostly concentrated in oasis areas with frequent human activities, while the highest habitat 
quality was generally concentrated in mountains and water bodies. 


3.2 Heterogeneity of habitat quality in the TRB 


Based on the spatial distribution of habitat quality of the TRB in the fourth periods (1990, 2000, 
2010 and 2018) and according to the calculation process of Equation 4, we tested the spatial 
autocorrelation of habitat quality using GeoDa software. The results shown in Figure 4 
demonstrated that when P<0.05, the global Moran's J index values of habitat quality during the 
four periods from 1990 to 2018 were 0.8465, 0.8437, 0.8429 and 0.8395, respectively, thereby 
confirming the existence of spatial autocorrelation of habitat quality, i.e., signs of heterogeneity. 
Generally, the global autocorrelation index displayed a continuous weakening trend between 1990 
and 2018, but the decreasing range was very small, which indicated that the spatial heterogeneity 
had a weak decreasing trend. 

We divided the spatial heterogeneity units of habitat quality in the TRB into three types 
according to Equation 5: cold spot areas, hot spot areas and random areas. Habitat quality of the 
TRB from 1990 to 2018 was then analyzed using the Getis-Ord G; index in ArcGIS; the results 
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of which are shown in Figure 5. According to ArcGIS statistics, cold spot areas accounted for 
42.03%, 41.67%, 54.27% and 63.97% of the study area in 1990, 2000, 2010 and 2018, 
respectively; hot spot areas accounted for 26.07%, 25.79%, 22.89% and 20.25%, respectively; 
and random areas accounted for 31.90%, 32.54%, 22.84% and 15.78%, respectively. Spatially, the 
distribution pattern of cold spot areas and hot spot areas of habitat quality was "cold inside and 
hot outside, cold in the east and hot in the west" and displayed obvious heterogeneity. Cold spot 
areas were mainly concentrated in the Taklimakan Desert and Kumtag Desert, and the KRB and 
CRB had more cold spot areas than the other sub-basins. The land use types in cold spot areas 
were mainly alpine desert, bare rock land or grassland with low vegetation coverage; therefore, 
habitat quality was low, with average habitat quality of 0.1240. The hot spots were mainly located 
in Oasis areas, mountain areas and mainstream of the Tarim River (ecological protection area). 
Due to the relatively high vegetation coverage in these areas, there was generally a high habitat 
quality (average habitat quality of 0.6594). Regarding time change, the range of cold spot areas 
exhibited a continuously increasing trend from 1990 to 2018, whereas hot spots in mountainous 
areas showed an increasing trend. 


3.3 Analysis of the change in Hemeroby index 


Based on the spatial distribution of land use types in the Tarim River Basin between 1990 and 
2018, we calculated the spatial distribution of Hemeroby index in the TRB using Equation 4; 
results are shown in Figure 6. According to the statistics (Table 6), Hemeroby index of the TRB 
during the four periods (1990, 2000, 2010 and 2018) were at high levels: 0.6459, 0.6472, 0.6586 
and 0.6608, respectively. This showed a general upward trend, which indicated that the impact of 
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Fig.4 Global Moran's / index statistics of habitat quality in the TRB in 1990 (a), 2000 (b), 2010 (c) and 2018 (d) 
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Fig. 5 Spatial distribution of cold spots and hot spots of habitat quality in the TRB in 1990 (a), 2000 (b), 2010 
(c) and 2018 (d) 
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Fig. 6 Distribution of Hemeroby index degree levels in the TRB in 1990 (a), 2000 (b), 2010 (c) and 2018 (d) 
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Table 6 Hemeroby index at different altitudinal gradients in 1990, 2000, 2010 and 2018 


Hemeroby index 


Altitudinal gradient 


1990 2000 2010 2018 
Gradient 1 0.6791 0.6810 0.6906 0.6905 
Gradient 2 0.6337 0.6331 0.6439 0.6486 
Gradient 3 0.5519 0.5522 0.5772 0.5820 
Gradient 4 0.6905 0.6486 0.5820 0.5950 


TRB 0.6459 0.6472 0.6586 0.6608 
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human disturbance on the ecological environment was intensifying between 1990 and 2018. In 
terms of spatial distribution, Hemeroby index was similar to the distribution pattern of cold spots 
and hot spots, displaying a pattern of "high inside and low outside, high in the east and low in the 
west". In addition to the desert areas, Hemeroby index with high level was mainly distributed in 
low-altitude areas (gradient 1 and gradient 2). These areas were mainly oases with frequent 
human activities; therefore, Hemeroby index was high. Hemeroby index with middle level was 
mainly distributed in the mountainous areas (gradient 3) with less human activities, more 
precipitation and higher vegetation coverage. There were few areas with low level of Hemeroby 
index, and they were mainly distributed around lakes and other water bodies. In terms of time 
scale, the areas of low- and middle-level ecological interference continued to decrease (with the 
exception of desert areas) between 1990 and 2018, and the ecological interference degree of the 
same altitude area displayed an increasing trend, showing that the impact of human activities on 
the ecological environment was intensifying. 


3.4 Linear relationship between habitat quality and disturbance factors 


3.4.1 Selection of disturbance factors 


Results in Section 3.3 showed that from 1990 to 2018, Hemeroby index in the TRB displayed an 
increasing trend, and habitat quality was in a state of degradation. As the premise and foundation 
of ecosystem service functions, habitat quality has a number of disturbance factors (Fellman et 
al., 2015; Hillard et al., 2015; Han et al., 2020). Hemeroby index is a comprehensive indicator of 
human disturbance and reflects the disturbance degree of landscape type (or land use type) change 
on an ecological environment. Therefore, during the analysis of habitat quality and disturbance 
factors, it is also necessary to consider other natural and human factors, in addition to the 
availability of data. Relevant studies (Sun et al., 2019; Aneseyee et al., 2020) have shown that 
population density, socioeconomic factors, topography and geomorphology have impacts on the 
distribution of habitat quality. The water system structure of the TRB is separate, and "water is 
oasis, no water is desert" is the typical characteristic of the basin. Water plays a crucial role in the 
ecological environment of the basin (particularly the mainstream of the Tarim River). Based on 
the above characteristics, this study chose Hemeroby index, population density, GDP, altitude and 
distance from water source (DWS) as the main disturbance factors. The above five variables were 
diagnosed using collinearity (Fig. 7), and the results demonstrated that the variance inflation 
factor (VIF) values for each disturbance factor in different years were less than 5.0, which 
indicated that there was no collinearity between the factors. 
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Fig. 7 Results of variance inflation factor (VIF) in collinearity diagnostics for the main disturbance factors. 
Cold, cold spots; Hot, hot spots; Random, random areas. H, Hemeroby index; PD, population density; GDP, gross 
domestic product; DWS, distance from water source. 
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3.4.2 Linear relationship between habitat quality and disturbance factors 

The correlation between habitat quality and disturbance factors in the three heterogeneous units 
can be seen in Figure 8. Different disturbance factors had different correlations with habitat 
quality in different years. A significant and negative linear relationship was found between habitat 
quality and Hemeroby index in the three heterogeneous units (P<0.01), and the correlation 
coefficient was largest, which indicated that Hemeroby index was the key disturbance factor that 
caused the degradation of habitat quality. Habitat quality showed a positive linear relationship 
with population density and GDP in cold spot areas and random areas and a significant and 
negative linear relationship in hot spot areas, which indicated that population density and GDP 
facilitated habitat quality degradation in hot spot areas but not in cold spot areas and random 
areas. The reason for this difference is that hot spot areas were the main population gathering 
regions with frequent human activities, which can directly or indirectly affect the habitat. 


(a) Cold spots (b) Hot spots (c) Random areas 


DWS + 


Altitude 


GDP 


PD 


-1.0 -0.5 0.0 0.5 -1.0 -0.5 0.0 0.5 -1.0 -0.5 0.0 0.5 
Correlation coefficient 


Fig. 8 Relationship between habitat quality and disturbance factors in cold spot areas, hot spot areas and random 
areas. *, statistically significant at the 5% level; **, statistically significant at the 1% level. 


There was a significant and positive correlation between habitat quality and altitude, indicating 
that habitat quality improved with increasing altitude, but the correlation coefficient in hot spot 
areas decreased between 1990 and 2018. There was a significant and negative linear relationship 
of between habitat quality and DWS in cold spot areas and random areas, and the correlation 
coefficient increased (from 1990 to 2018), indicating that habitat quality worsened in these areas 
farther away from water sources. Further, there was a significant and negative linear relationship 
between habitat quality and DWS in the period of 1990-2000 and then a significant and positive 
correlation in the period of 2010-2018, which indicated that the ecological environment of the 
TRB improved following the implementation of ecological management projects (particularly the 
ecological water conveyance project) after 2000. 

3.4.3 Spatial relationship between habitat quality and disturbance factors 

The spatial relationship between habitat quality and disturbance factors (Hemeroby index, 
population density, GDP, altitude and DWS) of the three heterogeneous units (cold spot areas, hot 
spot areas and random areas) in the TRB between 1990 and 2018 was statistically tested and 
compared using the regression modeling tool of the GeoDa software to enable the selection of the 
optimal spatial regression model. The evaluation indices of the spatial regression model included 
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correlation coefficient (R*), log (Likelihood) (logL), Akaike information criterion (AIC) and 
Schwartz criterion (SC). The range of R? values was 0.0000—1.0000, and the closer the R? values 
were to 1.0000, the better the regression effect of the spatial regression model. In addition, the 
larger the log value and the smaller the AIC and SC values were, the better the regression effect 
of the model. According to the OLS, it was necessary to judge the significance of the Lagrange 
multiplier (LM) and Robust lagrange multiplier (RLM). The larger the values of the two were, the 
better the regression effect of the spatial regression model. Table 7 shows that the global Moran’s 
I (error) values of the three heterogeneous units between 1990 and 2018 were significant at 
P<0.01 level, which indicated that there was a spatial dependence in all spatial regressions. 
Furthermore, LM (lag) and LM (error) all passed the significance test of 1%. Therefore, SEM or 
SLM regression was better suited to the study of habitat quality and disturbance factors than the 
OLS regression. The regression model with the best fitting effect was chosen following a 
comparison between the LM and RLM values of the three heterogeneous units for each year. By 
examining the advantages and disadvantages, the SLM was deemed to be more suitable for the 
spatial regression analysis. 


Table 7 Comparison of statistical test goodness of the spatial regression models in different heterogeneous units 
in 1990, 2000, 2010 and 2018 


Heterogeneous unit Statistic 1990 2000 2010 2018 
R 0.9343 0.9259 0.5232 0.4690 
logL 198,929.0 198,940.0 59,972.3 50,485.7 
AIC —397,845.0 —397,868.0 —119,933.0 —100,959.0 
SC —397,791.0 —397,813.0 —119,877.0 —100,904.0 
Cold spot areas Moran’s J (error) 68.9%" 72.0” 347.3" 388.0" 
LM (lag) 556.4" 415.0" 84,837.9"° 116,775.9" 
RLM (lag) 1.2 51.6" 94.3" 1128.1" 
LM (error) 4736.6" 5182.9" 120,561.3° 150,482.0° 
RLM (error) 4181.4" 4819.5" 35,817.7" 34,834.2"° 
R 0.8447 0.8367 0.3661 0.3367 
logL 47,430.2 45,911.6 17,201.2 15,911.7 
AIC —94,848.3 —91,811.2 —34,390.4 —31,811.3 
SC —94,796.2 —91,759.1 —34,338.4 —31,759.4 
Hot spot areas Moran’s 7 (error) 245.0" 242.3" 192.2*° 188.8” 
LM (lag) 3771.4" 4116.7" 20,683.1°° 21,503.1™ 
RLM (lag) 3638.0" 3501.3" 1402.6" 778.3" 
LM (error) 59,975.4"° 58,627.1°° 36,908.3™" 35,592.3™ 
RLM (error) 59,842.0°° 58,011.8" 17,627.7" 14,867.5"° 
R? 0.8928 0.8826 0.4390 0.3832 
logL 73,152.8 72,003.3 30,106.2 13,218.4 
AIC —146,294.0 —143,995.0 —40,200.4 —26,424.8 
SC —146,239.0 —143,940.0 —40,147.8 —26,374.4 
Random areas Moran’s J (error) 312.3" 319.7" 220.1" 175.3" 
LM (lag) 23,710.9"° 26,155.5"° 30,901.5" 20,351.3™ 
RLM (lag) 635.5" 615.1” 219.4" 148.4" 
LM (error) 97,456.1"" 102,120.4™ 48,391.0" 30,678.2"° 
RLM (error) 74,380.7"° 76,580.0°"" 17,709.0"" 10,475.3"° 


Note: R?, Correlation coefficient; logL, log(Likelihood); AIC, Akake information criterion; SC, Schwartz criterion; LM, Lagrange 
multiplier; RLM, Robust lagrange multiplier; **, statistically significant at the 1% level. 


According to the spatial regression coefficient (Table 8), the regression coefficients of habitat 
quality with Hemeroby index, altitude and DWS in each heterogeneous unit between 1990 and 
2018 was negative (P<0.01), which indicated that there was a significant negative correlation 
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between habitat quality and these disturbance factors. However, the correlation was different in 
the three heterogeneous units. The negative correlation coefficient between habitat quality and 
Hemeroby index was the largest, which indicated that Hemeroby index had the greatest impact on 
habitat quality. The impact of Hemeroby index on habitat quality in cold spot areas and random 
areas was greater than that in hot spot areas. There was a significant and positive correlation 
between habitat quality and GDP (P<0.01), and the correlation was greater in hot spot areas. 
There was a positive correlation between habitat quality and population density in cold spot areas 
and random areas, but there was a negative correlation between them in hot spot areas, indicating 
that there was a greater likelihood that human aggregation in hot spot areas will affect habitat 
quality. Generally, the effects of habitat quality in different heterogeneous units on the spatial 
dependence of disturbance factors were different. 


Table 8 Spatial regression results of habitat quality with disturbance factors in different heterogeneous units in 
1990, 2000, 2010 and 2018 based on the spatial error model 


Heterogeneous unit Variable 1990 2000 2010 2018 
Constant 1.8829" 1.8855” 1.2565" 1.1626" 
H —2.4739"* —2.4775™ —1.5976" —1.4406™ 
PD 0.0013 0.0051” 0.0014” 0.0043” 
Cold spot areas . ; , 
GDP 0.0068% 0.0349% 0.0571” 0.0158" 
Altitude —0.0068"™ —0.0069" —0.0016™ —0.0017™ 
DWS —0.0073™ —0.0075™ —0.0019™ —0.0041™ 
Constant 1.6167" 1.6592" 1.3339" 1.3222™ 
H —1.901™ —1.8900"* —1.2942™ —1.2416" 
PD —0.0028"* —0.0086 —0.0045™ —0.0057"* 
Hot spot areas Pe . . ‘a 
GDP 0.0169" 0.0580" 0.0876” 0.0303” 
Altitude —0.0025™ —0.0026™ —0.0017™ —0.0021™ 
DWS —0.0019™ —0.0021™ —0.0034™ —0.0043™ 
Constant 1.8317" 1.8260" 1.2711" 1.2653" 
H —2.3133™ —2.2993" —1.4637" —1.3628" 
PD 0.0117" 0.0097" 0.0048" 0.0071” 
Random areas ~ oe Pv ov 
GDP 0.0146 0.0652” 0.0935" 0.0103" 
Altitude —0.0016™ —0.0016™ —0.0047™ —0.0016™ 
DWS —0.0023™ —0.0025™ —0.0074™ —0.0010™ 


Note: H, Hemeroby index; PD, population density; GDP, gross domestic product; DWS, distance from water source. 
significant at the 1% level. 


, Statistically 


4 Discussion 


4.1 Spatial and temporal characteristics of habitat quality heterogeneity 


This study discovered that habitat quality of the TRB displayed a continuous degradation trend 
between 1990 and 2018, and the sub-basins (with the exception of the HRB) also showed 
different degrees of degradation. There were differences in habitat quality in different sub-basins 
and at different altitudes. Habitat quality was also generally better in the northern region of the 
study area than in the southern region and in high-altitude areas than in low-altitude areas. This 
could be because there was a greater vegetation coverage in the northern region than in the 
southern region and the frequency of human activities in high-altitude areas was lower. These 
results were consistent with those obtained by Liu et al. (2020) on the spatial and temporal 
changes in habitat quality in Xinjiang of China: the overall habitat quality in Xinjiang displayed a 
degradation trend, and the areas with high values of habitat quality were mainly located at the 
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edge of mountains and basins (Liu et al., 2020). At the same time, the results of our study were 
also similar to those obtained by research on the landscape biodiversity of the YRB in the TRB by 
Huang et al. (2020): areas with low values of habitat quality were mainly distributed in regions 
with a poor natural environment and areas with high values of habitat quality were mainly 
distributed in natural forests and grassland reserves (Huang et al., 2020). Our study also 
confirmed that there was a significant spatial heterogeneity in habitat quality of the TRB, and the 
study area can be divided into three heterogeneous units in terms of habitat quality: cold spot 
areas, hot spot areas and random areas (Fig. 5). In the three heterogeneous units, cold spot areas 
were mainly concentrated in low-altitude areas, which were the main regions of human activities, 
in addition to the two large deserts (Taklamakan Desert and Kumtag Desert). Hot spot areas were 
mainly distributed in high-altitude mountainous regions or ecological protection regions, as these 
areas had relatively high vegetation coverage. 


4.2 Effects of different disturbance factors on habitat quality 


Hemeroby index quantifies the disturbance degree of human activities on an ecosystem by 
considering the degree of disturbed landscape or habitat in relation to the natural landscape or 
habitat (Anselin, 1988; Chen et al., 2010). Cultivated land, construction land and unused land 
were threat source factors to habitat quality. Their changes significantly impacted habitat quality. 
Figure 9 showed that the proportion of threat source factors in the three heterogeneous units was 
largest (with the exception of natural vegetation), which explained the significant and negative 
correlation between habitat quality and Hemeroby index. However, the negative correlation 
between habitat quality and Hemeroby index has decreased since 2000, potentially because the 
comprehensive management policy of the TRB has effectively restrained the deterioration of the 
ecological environment in the region. Among other disturbance factors, population density and 
GDP exhibited opposite correlations in cold spot areas, random areas and hot spot areas. 
Combined with the gradual deterioration of habitat quality in the oasis areas in Figure 3 and the 
increasing proportion of cultivated land in Figure 9, this showed that human activities had a 
negative impact on habitat quality of hot spot areas. 


o 
= 
RS Cold spots Hot spots Random areas a8 
v I 
S 
9 
o 
Area percentage (%) Area percentage (%) Area percentage (%) 


Cultivate land gg 9.0000} 0.0000} 0.1698 3.3306} 3.4756} 6.9095 
Forestland M 0.0021] 0.0042} 0.0857) 3.9638] 4.0679} 4.1679 0.7045} 0.7878} 1.103 
Grasland M8 0.0104} 0.0042) 5.1956 75.2913] 74.2545 |74.5015 18.0500] 16.7512|14.3055 
Construction land $ 0.0000} 0.0000) 0.1828) 0.0000} 10.6822) 5.9701 1.7172] 1.8319} 2.01 
Water body gg 0.0021} 0.0021} 0.3284 0.1381} 0.1328} 0.1725 0.3825} 0.3804) 0.6653 
Unused land = 99.9854 |99.9895|94.0377 6.7017] 7.3870} 8.2784 74.3327|75.0175|73.7925 


4.8131} 5.2312) 8.1141 


Fig. 9 Area percentages of land use types in the three heterogeneous units. The arrow direction and years 
indicated that the circles from outside to inside were 1990, 2000, 2010 and 2018, respectively. Values in the table 
represented the area percentages of land use types. The colors for area percentages of land use types were 
consistent with those of years. 


4.3 Limitations of the nVEST model 


The InVEST model compensates for the shortcomings of traditional ecosystem service 
assessment methods and provides an effective way for assessing ecosystem services. In recent 
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years, as a function module of the InVEST model, habitat quality assessment has been widely 
conducted (Sharp et al., 2016; Babbar et al., 2021). However, the model has some shortcomings 
when assessing habitat quality. For example, it obtains the habitat quality in a study area by 
accumulating the influence of threat source factors on habitat quality, but the simple accumulation 
of threat source factors is not exactly equal to the comprehensive impact of these stress factors. 
Therefore, the model requires further improvement in terms of the usage process. The parameters 
of threat source factors and habitat sensitivity are subjective and require further study in the 
future. 


5 Conclusions 


Through an analysis of the spatial and temporal changes in habitat quality in the TRB between 
1990 and 2018, this study discussed the spatial heterogeneity of habitat quality and quantified its 
relationship with interference factors, thereby providing a new comprehensive method for the 
quantification of the impact of threat source factors on habitat quality. The main conclusions are 
as follows. 

Temporally, the overall level of habitat quality in the TRB was poor and in a state of continuous 
degradation. Spatially, habitat quality in the northern region was greater than that in the southern 
region and greater in high-altitude areas than in low-altitude areas. The results showed that habitat 
quality exhibited spatial heterogeneity in the TRB, represented by cold spot areas, hot spot areas 
and random areas. There was a significant negative linear relationship between habitat quality and 
Hemeroby index in the three heterogeneous units, and the correlation coefficient was greater than 
the correlation coefficients of habitat quality with population density, GDP, altitude and DWS. 
There was a spatial dependence between habitat quality and disturbance factors in the three 
heterogeneous units. During the four periods (1990, 2000, 2010 and 2018), habitat quality in each 
heterogeneous unit had a significant and negative correlation with Hemeroby index, altitude and 
DWS, but the correlation coefficient with Hemeroby index was the largest, which indicated that 
Hemeroby index had the greatest impact on habitat quality. 
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